Nonequilibrium dynamical ferromagnetism of interacting Single-Molecule Magnets 
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We propose a nonequilibrium Monte Carlo (MC) approach to explore nonequilibrium dynamical 
ferromagnetism of interacting single molecule magnets (SMMs). Both quantum spin tunneling and 
thermally activated spin reversal are successfully implemented in the same MC simulation frame- 
work. Applied to a typical example, this simulation method satisfactorily reproduces experimental 
magnetization curves with experimental parameters. Our results show that both quantum and clas- 
sical effects are essential to determine the hysteresis behaviors. This method is effective and reliable 
to gain deep insights into SMMs. 
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Introduction. The single-molecule magnet (SMM) is 
very interesting because of its potential applications 
in spintronic devices 0, 0], quantum computing in- 
formation storage optical transistors 0, and so on. 
The most famous examples are Mni2 0, 13] , Fes § , and 
Mn4 0]. A SMM usually has well-defined spin ground 
state and large magnetic anisotropy, and can exhibit hys- 
teresis loops at low enough temperatures 10] . Landau- 



Zener[lll [12| (LZ) model has been used to explain the 



steps of the hysteresis loops (lfj, [LJ, ll4|. On theoret- 
ical side, numerically exact method has been used to 
the phenomenological spin model for the Fes and Mni2 
systems, and calculated results show that LZ model is 
reliable for usual experimental field sweeping rates in 
describing the two-level problems at the level crossing 
fields 



15 



161 ] , but experimental hysteresis loops cannot be 
explained because thermally activated processes are not 
considered [Io|. Thermally activated processes, however, 
play important roles in determining the nonequilibrium 
magnetization dynamics and even level populations. A 
quantitative method which is able to reproduce the full 
hysteretic magnetization dynamics of SMMs is in need. 

Here, we propose a nonequilibrium Monte Carlo 
(NQMC) approach taking both the LZ quantum tunnel- 
ing effect and the thermally activated processes into ac- 
count, and hereby satisfactorily reproduce and explain 
experimental magnetization behaviors of SMMs. 

Spin model. For simplicity without losing physics, we 
arrange SMMs on a finite two-dimensional (2D) rectan- 
gular lattice of Li x L2 , and consider inter-SMM magnetic 
dipolar and exchange interactions. Actually, some single 
layers of patterned SMMs has already been grown on 
specific solid sur 

faces USQUE!, 

and it has been proved 
that a SMM, when put on a Au surface, can still keep its 
essential properties [4j. We use giant spin approximation. 
Every SMM has a spin of 5. We use a typical lattice 
constant a = 1.5 nm, which is an intermediate value of 
well-known SMMs. The model Hamiltonian reads 
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s, 



(1) 



where jfj and Jij describe the magnetic dipolar and ex- 
change interactions between spins at i and j. The factor 
1/2 is due to the double counting in the summation. The 
first term describes the part for all the single SMMs, and 
Hf is given by 



D 2 (Sfy - D^Slf + Hf + 9HBB.SI , (2) 



where D2 and D4 are positive anisotropic parameters, g 
the Lande g- factor (here g = 2 is used), \ib the Bohr mag- 
neton, B z (t) the external magnetic field in the z direc- 
tion. Si = (Sf,Sf, Sf) is the quantum spin operator for 
the ith SMM. As for the transverse term H\ v , it is usually 
defined as Hf = E[{Sff - (Sf) 2 } + g/j, B B x Sf, where E 
is the second-order transverse anisotropic constant and 
B x is the transverse external field in the x direction. As 
usual, B x is a constant and B z (t) is a linear function of 
the time t with the sweeping rate v. 

LZ tunneling probability. The easy axis determines two 
equilibrium spin directions in the iz direction. When 
calculating the LZ tunneling probability, we treat the 
SMM spins coupled with the ith spin by using the mean 
field approximation. As a result, for the zth SMM, we 
have the effective one-SMM Hamiltonian 



H i = H2+gn B BfSf 
where the effective mean field B° s is given by 
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eff 
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Ji 



(3) 



(4) 



where S e ^ is the equilibrium value [S or —S) for the ith 
spin along the easy axis. This approximation is natu- 
ral in the NQMC simulation, as will be clarified in the 
following. Then the standard diagonalization technique 
can be used to solve Eq. ©. Hamiltonian Eq. has 
25+1 energy levels, which can be labeled by the quantum 
numbers m = S, S—l, • • • , —(5—1), —5 to a first-order ap- 
proximation. If without the transverse part Hf r , Eq. ([3]) 
is diagonal, and there are level crossings at some special 
field values. When the transverse term H*" 1 is taken into 
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account, the level crossings become avoided [lol, 13, 14 1. 
When B z is swept close to B m<m i at which the avoided 
level crossing happens between states to and m', quan- 
tum tunneling occurs between the two states. Strictly 
speaking, this tunnel ing is beyond LZ model which is a 
two-states theory [HI [l2j. but LZ tunneling can be used 
as a good approximation for each of the tunneling pro- 
cesses because LZ transition time is usually very short 
compared to time increments between the two successive 
avoided level crossings HI, 3, 2(j • The nonadiabatic LZ 
tunneling probability P m ,m' is given by[ll|, [T^ 



exp 



7rA 



2 



2hgfiB\m 



(5) 



where the tunnel splitting A mjm > describes the energy 
gap at the avoided crossing of states to and to'. B mrn i 
and A m , m > can be calculated by diagonalizing Eq. ([3]). 

Thermally activated reversal rate. In order to calculate 
thermal spin reversal rate, we use classical spin approxi- 
mation which is reasonable and reliable because the spin 
usually is very large for SMMs. As a result, the energy 
of the ith SMM is given by 



Ei = -D 2 (S:, 



z\2 



z\4 



htSf 



(6) 



where hi = q[Ib(B z + Bf ). Each of the spins has two 
equilibrium orientations (±5) along the easy axis. We 
use the angle Oi to describe the ith spin's deviation from 
its original equilibrium orientation. All the other an- 
gle values (0 < 6i < it) are the transition states for 
the spin to reverse its orientation. We express Sf as 
5* q cos Oi . Usually, there exists a maximum in the curve 
of Ei(cos9i), and the maximum determines the energy 
barrier for the spin reversal 21, 22| . Defining x\ = cos Oi, 
we have — 1 < Xi < 1. Then the barrier is given by 



Ei(xi 
AEi = { Ei ( Xi 
Ei(xi 



x?) 

-1) 
1) = 







if \x^\ < 1 
if x? < -1 
if xT > 1 



where x™ (according to the maximum) is defined as 



(7) 



(8) 



where A d = (q/2) 2 + (p/3) 3 , p = D 2 /(2D 4 S 2 ), and 
q = -/i l 5 l cq /(4D 4 5 4 ). Then the reversal rate of the 
ith spin is determined by the Arrhenius law 23], Ri — 
Rq exp(— AEi/ (ksT)), where ks is the Boltzmann con- 
stant and Rq the characteristic attempt frequency. 

Simulation method and parameters. To set up our 
NQMC simulation, we define MC steps by the time 
points, t n = At ■ n, where n nonnegative integers. Every 
spin takes cither S or —S at each of the times t n . A spin 
can be reversed by overcoming the thermal barrier and 
tunneling through LZ mechanism. The probability P, for 




FIG. 1: (color online). Simulated magnetization curves 
{M/Ms vs. B z in T) with three sweeping rates: v = 0.002 
(a), 0.02 (b), and 0.14 (c) T/s. For all the three panels, the 
temperatures are T = 0.10 (•), 0.20 (■), 0.30 (a), 0.40 (T), 
0.50 (♦), 0.60 (X), 0.70 (O), and 1.00 (□) K. 



a thermally activated spin reversal within a MC step is 
defined by P l = At 



Ri 21 



22] 
©■ 



and the LZ tunneling 
Actually, we have the 



probability is given by Eq. 
magnetization m = S at the beginning of field sweeping, 
and then need to consider only the avoided crossings at 
(5, to'), where ml = —S, —S +!,-■■. When the field is 
swept to Bs t m'i we consider the LZ tunneling by adding 
the LZ probability Ps,m> to the thermally activated re- 
versal probability. Because our Monte Carlo time step is 
much larger than the LZ transition and spin relaxation 
times, a spin will have long enough time to transit to the 
possible lowest level at current field after a LZ tunnel- 
ing is finished 24|. Therefore, our NQMC simulation is 
self-consistent. 

For convenience in comparing with experiment, we 
choose the Mn 4 system to demonstrate the NQMC 
method and its power. Hence we have S = 9/2, Z^/fes = 
0.608K, D 4 /k B = 3.8mK, and E/k B = 32mK from Ref. 
(3- We assume that our 2D lattice is put on a metal 
surface and hence use a carrier-mediated 2D interaction 
for Jij (J/ks — 5.8mK for the nearest spins) 25]. This 
inter-SMM interaction, weak enough, does not destroy 
the essential SMM properties. In addition, we already 
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show that simulated results are not sensitive to specific 
choices of the interaction, and other interactions yield 
similar results. We take B x = 0.02T, At = 1ms, and 
i? = 10 9 /s. Such parameter choices guarantee the good 
balance between computational demand and precision. 
In our simulations, field B z is swept from -2.5T to 2.5T 
in the forward process, and the full hysteresis loop is 
obtained simply by using the loop symmetry. Every hys- 
teresis loop is calculated by averaging over 100 runs to 
reduce statistical errors. 

Simulated results and discussions. We have done 
our NQMC simulations with various {L\,L<i) values and 
experimentally accessible sweeping rates v (0.002—0.14 
T/s) and temperatures T (0.10-1.00 K)[H]. Presented 
in Fig. 1 are our simulated magnetization M (normalized 
to the saturated value M s ) curves against the applied 
sweeping field B z with L\ — L 2 = 10 for three sweep- 
ing rates, namely 0.002, 0.02, and 0.14 T/s. For each 
of the three cases, we present the results for eight tem- 
peratures: 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, and 
1.00 K. It is clear that there exist hysteresis loops for low 
enough temperatures and complete paramagnetic mag- 
netizations are obtained for high enough temperatures. 
We can clearly see four steps between and 2T for low 
enough temperatures. In addition, we have done NQMC 
simulations for various values of L\ and L 2 . For square 
lattices (L = L\ = L 2 ), the step heights decrease a little 
with increasing L from 8 to 50. For rectangular lattices 
with L\ x L 2 = 24 2 , the step heights change a little for 
different L\ and L 2 . It is clear that the step structures 
are far from those predicted by LZ tunneling alone. 

In fact, thermally activated reversals play important 
roles in determining the magnetization curves. Our sim- 
ulations show that the number of the steps decreases with 
increasing the temperature for given value of the sweep- 
ing rate v. It can be seen that the area enclosed by a 
hysteresis loop increases monotonically with increasing 
the sweeping rate and decreasing the temperature, and 
the smallest visible loop is at 0.70 K for v = 0.14 T/s, 
0.60 K for 0.02 T/s, and 0.50 K for 0.002 T/s, respec- 
tively. The thermal effect increases substantially when 
we increase the temperature from 0.3K on. There will 
be no step structures for usual sweeping rates when the 
temperature is higher than 0.5K. The hysteresis loop will 
disappear when the temperature is too high. 

Conclusion. We have proposed a NQMC approach to 
explore the nonequilibrium dynamical ferromagnetism of 
interacting SMMs. Both quantum spin tunneling and 
thermally activated spin reversal are successfully imple- 
mented in the same Monte Carlo simulation framework. 
Applied to the Mn4 system, this simulation method sat- 
isfactorily reproduces experimental magnetization curves 
with experimental parameters. Our results show that 
both quantum and classical effects are essential to de- 
termine the hysteresis loops in the magnetization curves. 
This method is effective and reliable, and can be used to 



gain deep insights into SMMs and other quantum nano- 
magnets with strong anisotropy. 
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